Recocimiento simulado¶

Supongamos que queremos minimizar una función $f:\mathbb{R}\rightarrow \mathbb{R}$ y tenemos un kernel de transición $Q$ para realizar una búsqueda aleatoria de este mínimo.

Algoritmo de busqueda aleatoria cruda¶

1) Sea $x_0$ un candidato inicial del mínimo de $f$ y $y_0=f(x_0)$

2) Para $i$ desde $1$ hasta $N$

2.1) Simule $x_n\sim Q(\cdot|x_{n-1})$ y tome $y_n=f(x_n)$

3) Tome $y^\star=\min\left\{y_0,\ldots , y_N\right\}$ como aproximación al valor que minimiza $f$.

El algoritmo de recocimiento simulado genera una cadena en la cual se acepta la propuesta cuando $f(x_{\text{prop}}) \leq f(x_{\text{prev}}) $ En caso de que $f(x_{\text{prop}}) > f(x_{\text{prev}})$ se considera un paso de aceptación tipo Metropolis-Hastings haciendo uso de la función de densidad de probabilidad

$$ p(x)=\frac{e^{-\frac{f(x)}{T}}}{c_p} $$

donde $c_p$ es la constante de normalización y el parámetro $T$ tiene la interpretación de temperatura en el contexto de física estadística; y una probabilidad de aceptación $$ \alpha(x_{\text{prop}}|x_{\text{prev}})=\frac{ p(x_{\text{prop}}) }{ p(x_{\text{prev}}) }= e^{-\frac{\left( f(x_{\text{prop}})-f(x_{\text{prev}})\right)}{T}} \mathbb{1}_{\left\{ f(x_{\text{prop}}) > f(x_{\text{prev}}) \right\}}. $$

Por lo que con probabilidad positiva en el caso de que $ f(x_{\text{prop}}) > f(x_{\text{prev}}$ podemos aceptar $x_{\text{prop}}$ a pesar de que este valor es un peor candidato para la minimización de $f$. Observemos que si la temperatura $T$ es cercana a cero entnces la probabilidad de aceptación $\alpha$ será cercana a cero. Esto puede ser de utilidad para salir de regiones cercanas a un mínimo local cuando la temperatura $T$ no es cercana a cero.

Algoritmo de recocimiento simulado¶

1) Sea $x_0$ un candidato inicial del mínimo de $f$ y $y_0=f(x_0)$

2) Para $i$ desde $1$ hasta $N$

2.1) Simule $x_{\text{prop}}\sim Q(\cdot|x_{n-1})$ y tome $y_{\text{prop}}=f(x_n)$

2.2) Si $y_{\text{prop}} \leq y_{i-1}$ tome $x_i=x_{\text{prop}}$ y $y_i=y_{\text{prop}}$. Proceda a 2.3) sólo en el caso contrario $y_{\text{prop}} > y_{i-1}$

2.3) Tome $U\sim \text{Uniforme}(0,1)$

2.3.1) Si $U < \alpha(x_{\text{prop}}|x_{i-1})$ entonces tome $x_i=x_{\text{prop}}$ y $y_i=y_{\text{prop}}$; en caso contario, $U\geq \alpha(x_{\text{prop}}|x_{i-1})$ tome $x_i=x_{i-1}$ y $y_i=y_{i-1}$

A lo largo del algoritmo anterior es recomendable calcular el mínimo de $f$ sobre la cadena $x_0,x_1,\ldots , x_N$ construida.

In [131]:
# Función para propuesta con cambio por entrada
function neighbor(x₀::Array{Float64,1},i::Int64,v::Array{Float64,1})
    n = length(x₀)
    return (x₀ + v[i]*rand(Uniform(-1,1))*I(n)[i,:]).nzval
end
Out[131]:
neighbor (generic function with 1 method)
In [132]:
# Función para temperatura en la probabilidad de aceptación
temperature(t::Real) = 1 / log(t)
Out[132]:
temperature (generic function with 1 method)
In [133]:
# Algoritmo de recocimiento simulado con cambio por entrada
function recocimiento_simulado_entry(f::Function,t::Function,neigh::Function,x₀::Array{Float64,1},N::Int64,
            v::Array{Float64,1}=[0.0])
    n = length(x₀)
    if v ==  [0.0]
        v = ones(n)
    end
    x = x₀
    x_v = Array{Float64,1}[]
    push!(x_v,x)
    x_min = x₀
    y = f(x)
    y_min = y
    for j in 1:N
        t = temperature(j)
        for i in 1:n
            x_prop = neigh(x,i,v)
            y_prop = f(x_prop)
            if  y_prop <= y
            x = x_prop
            y = y_prop
            else
                p = exp( (y-y_prop)/t )
                if rand(Uniform()) < p
                    x = x_prop
                    y = y_prop
                end
            end
            
            if  y_prop <= y
                x = x_prop
                push!(x_v,x)
                if y_prop <= y_min
                    x_min = x
                end
            else 
                p = exp( (y-y_prop)/t )
                if rand(Uniform()) < p
                    x = x_prop
                    push!(x_v,x)
                else
                    push!(x_v,x)
                end
            end      
        end
        if y < y_min
            x_min = x
            y_min = y
        end
        push!(x_v,x)
    end
    return x_v[end], x_min
end
Out[133]:
recocimiento_simulado_entry (generic function with 2 methods)
In [140]:
# Función para propuesta
function neighbor(z)
  [rand(Uniform(z[1] - 1, z[1] + 1)), rand(Uniform(z[2] - 1, z[2] + 1))]
end
Out[140]:
neighbor (generic function with 2 methods)
In [134]:
# Algoritmo recocimiento simulado
function recocimiento_simulado(f::Function,t::Function,neigh::Function,x₀::Array{Float64,1},N::Int64)
    x = x₀
    x_v = Array{Float64,1}[]
    push!(x_v,x)
    x_min = x₀
    y = f(x)
    y_min = y
    for j in 1:N
        T = t(j)
        x_prop = neigh(x)
        y_prop = f(x_prop)
        if  y_prop <= y
            x = x_prop
            y = y_prop
        else
            p = exp( (y-y_prop)/T )
            if rand(Uniform()) < p
                x = x_prop
                y = y_prop
            end
        end
        if y < y_min
            x_min = x
            y_min = y
        end
        push!(x_v,x)
    end
    return x_min, x_v
end
Out[134]:
recocimiento_simulado (generic function with 1 method)
In [135]:
# Función de Rosenbrock con un mínimo en (1,1)
function rosenbrock(x::Array{Float64,1})
  (1 - x[1])^2 + 100(x[2] - x[1]^2)^2
end
Out[135]:
rosenbrock (generic function with 1 method)
In [136]:
# Implementación reocimiento simulado para Rosenbrock
recocimiento_simulado_entry(rosenbrock,temperature,neighbor,[0.1,0.1],10000)[1]
Out[136]:
2-element Array{Float64,1}:
 0.9867330094947029
 0.9629929384896192
In [144]:
gif(anim_simul_anneal, "anim_simul_anneal.gif", fps = 50)

┌ Info: Saved animation to 
│   fn = /home/workstation/Documents/Simulacion/anim_simul_anneal.gif
└ @ Plots /home/workstation/.julia/packages/Plots/lzHOt/src/animation.jl:104
Out[144]:
In [146]:
gif(anim_simul_anneal, "anim_simul_anneal.gif", fps = 50)

┌ Info: Saved animation to 
│   fn = /home/workstation/Documents/Simulacion/anim_simul_anneal.gif
└ @ Plots /home/workstation/.julia/packages/Plots/lzHOt/src/animation.jl:104
Out[146]:

Algoritmos genéticos¶

Los algoritmos genéticos son algoritmos de optimización vía busqueda aleatoria que hacen uso de una analogía con la manera en que los genes evolucionan. Para esto se consideran tres componentes:

1) Entrecruzamiento: Sea $\mathcal{C}$ el conjunto sobre el que deseamos optimizar. Dados dos puntos $c_1$, $c_2\in \mathcal{C}$, que son interpretados como cromosomas se tiene una función de entrecruzamiento, digamos $cross:\mathcal{C}\times \mathcal{C}\rightarrow \mathcal{C}$, que produce un nuevo punto que interpretamos como el "decendiente" de $c_1$ con $c_2$.

2) Mutación: Dado $c\in\mathcal{C}$ decimos que una mutación de $c$ está dada por $mut(c)$ con $ mut :\mathcal{C}\rightarrow \mathcal{C}$.

3) Aptitud: Dado $c\in\mathcal{C}$ decimos que tiene aptitud (fitness) $f(c)$ con $ f :\mathcal{C}\rightarrow \mathbb{R}$.

4) Selección: Dada una población $\pmb{c}\in\mathcal{C}^n$ podemos se seleccionar una subpoblación $s(\pmb{c})$ con $s:\mathcal{C}^n \rightarrow \mathcal{C}^m$, $m<n$, en la que usualmente se eligen con alta probabilidad los elementos de la población inicial más aptos.

In [16]:
# Ejemplo del vendedor viajero
function ciudades_sim( n, map_lim)
    c = Dict{String,Real}[]
    for c_run in 1:n
        push!(c, Dict(
                "id" => c_run, 
                "x" => round(rand() * map_lim), 
                "y" => round(rand() * map_lim) ) )
    end
    return c
end
# calling generate_cities method to create cities
ciud_dict = ciudades_sim(5,500)
Out[16]:
5-element Array{Dict{String,Real},1}:
 Dict("id" => 1,"x" => 228.0,"y" => 131.0)
 Dict("id" => 2,"x" => 94.0,"y" => 415.0)
 Dict("id" => 3,"x" => 127.0,"y" => 38.0)
 Dict("id" => 4,"x" => 278.0,"y" => 273.0)
 Dict("id" => 5,"x" => 458.0,"y" => 343.0)
In [48]:
function crossover( t1::Array{Int64,1}, t2::Array{Int64,1}, cross_i::Int64)
    tnew_1 = t1[1:cross_i]
    for c in tnew_1
        if c in t2
            c_loc = findfirst( i -> i == c, t2)
            splice!( t2, c_loc)
        end
    end
    return vcat( tnew_1, t2)
end
Out[48]:
crossover (generic function with 2 methods)
In [49]:
crossover( [1,2,3,4,5], [5,3,1,2,4], 2 )
Out[49]:
5-element Array{Int64,1}:
 1
 2
 5
 3
 4
In [50]:
function mutate(t::Array{Int64,1})
    mut_1 = rand(1:length(t))
    mut_2 = rand(1:length(t))
    t[mut_1], t[mut_2] = t[mut_2], t[mut_1]
    return t
end
Out[50]:
mutate (generic function with 2 methods)
In [51]:
mutate( [5,3,1,2,4] )
Out[51]:
5-element Array{Int64,1}:
 5
 3
 1
 2
 4
In [53]:
function fitness(c::Array{Int64,1},D::Array{Dict{String,Real},1})
    fit = 0
    c = vcat(1, c, 1)
    for loc in 1:length(c) - 1
        p1 = [ D[c[loc]]["x"], D[c[loc]]["y"] ]
        p2 = [ D[c[loc+1]]["x"], D[c[loc+1]]["y"] ]
        fit += euclidean(p1,p2)
    end
    return fit
end
Out[53]:
fitness (generic function with 2 methods)
In [84]:
fitness( [5,3,1,2,4], ciud_dict)
Out[84]:
1597.1842840869992
In [38]:
function inic_pop(n_pop::Int64,l_pop::Int64,D::Array{Dict{String,Real},1})
    c_vec = Dict{String,Any}[]
    for _ in 1:n_pop
        c = [1:l_pop;]
        for i in 1:length(c)
            u = rand(1:l_pop)
            c[i], c[u] = c[u], c[i]
        end
        push!(c_vec,  Dict( "id" => c, "fit" => fitness(c,D) ) )
    end
    return c_vec
end
Out[38]:
inic_pop (generic function with 1 method)
In [40]:
inic_c_pop = inic_pop(3,5,ciud_dict)
Out[40]:
3-element Array{Dict{String,Any},1}:
 Dict("id" => [5, 4, 2, 1, 3],"fit" => 1326.9704458779634)
 Dict("id" => [2, 3, 4, 5, 1],"fit" => 1477.7306971101777)
 Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
In [62]:
function evolve( inic_c::Array{Dict{String,Any},1}, n_ev::Int64, cross_i::Int64, D::Array{Dict{String,Real},1})
    for i in 1:n_ev
        n_pop = length(inic_c)
        c_dict = inic_c[i]
        u1 = rand(1:n_pop)
        u2 = rand(1:n_pop)
        t1 = copy( inic_c[u1]["id"])
        t2 = copy( inic_c[u2]["id"])
        new = crossover( t1, t2, cross_i)
        new = mutate(new)
        push!(inic_c, Dict( "id" => new, "fit" => fitness(new,D) ) )
    end
    return sort!(inic_c, by=x -> x["fit"], rev=false)
end
Out[62]:
evolve (generic function with 1 method)
In [65]:
evolve(inic_c_pop, 50, 2, ciud_dict)
Out[65]:
57-element Array{Dict{String,Any},1}:
 Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
 Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
 Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
 Dict("id" => [1, 4, 5, 2, 3],"fit" => 1230.4671534526028)
 Dict("id" => [3, 2, 4, 5, 1],"fit" => 1254.0912080872085)
 Dict("id" => [1, 5, 4, 2, 3],"fit" => 1254.0912080872085)
 Dict("id" => [3, 2, 4, 5, 1],"fit" => 1254.0912080872085)
 Dict("id" => [5, 4, 2, 3, 1],"fit" => 1254.0912080872085)
 Dict("id" => [3, 2, 4, 5, 1],"fit" => 1254.0912080872085)
 Dict("id" => [1, 3, 2, 4, 5],"fit" => 1254.0912080872085)
 Dict("id" => [5, 4, 2, 3, 1],"fit" => 1254.0912080872085)
 Dict("id" => [3, 4, 5, 2, 1],"fit" => 1294.8367593973094)
 Dict("id" => [2, 5, 4, 3, 1],"fit" => 1294.8367593973094)
 ⋮
 Dict("id" => [3, 4, 5, 1, 2],"fit" => 1550.6099349009326)
 Dict("id" => [3, 4, 5, 1, 2],"fit" => 1550.6099349009326)
 Dict("id" => [3, 4, 5, 1, 2],"fit" => 1550.6099349009326)
 Dict("id" => [4, 3, 5, 2, 1],"fit" => 1565.0505976063453)
 Dict("id" => [4, 3, 1, 2, 5],"fit" => 1565.0506101258716)
 Dict("id" => [2, 4, 3, 5, 1],"fit" => 1588.674652240951)
 Dict("id" => [5, 1, 2, 4, 3],"fit" => 1588.6746647604773)
 Dict("id" => [5, 1, 3, 4, 2],"fit" => 1588.6746647604773)
 Dict("id" => [3, 4, 2, 1, 5],"fit" => 1588.6746647604775)
 Dict("id" => [5, 3, 1, 2, 4],"fit" => 1597.1842840869992)
 Dict("id" => [2, 1, 4, 3, 5],"fit" => 1820.8237731099682)
 Dict("id" => [5, 3, 4, 1, 2],"fit" => 1820.8237731099684)
In [95]:
fit_run = Inf
c_run = Int64[]
for a in 1:5
    for b in deleteat!( [1:5;], a)
        for c in deleteat!( [1:5;], sort([a,b]))
            for d in deleteat!( [1:5;], sort([a,b,c]))
                for e in deleteat!( [1:5;], sort([a,b,c,d]))
                    fit_loc = fitness( [a,b,c,d,e], ciud_dict )
                    if fit_loc < fit_run
                        fit_run = fit_loc
                        c_run = [a,b,c,d,e]
                    end
                end
            end
        end
    end
end
In [96]:
c_run
Out[96]:
5-element Array{Int64,1}:
 1
 3
 2
 5
 4
In [97]:
fit_run
Out[97]:
1230.4671534526028